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1.  INTRODUCTION 

In  Reliability  studies,  one  of  the  basic  problems  is  the  estimation  of  the 
lifetime  distribution  of  a  product.  In  medical  research,  it  is  the  estimation  of  the 
lifetime  distribution  of  a  patient  under  remission  of  a  particular  disease.  For 
example,  in  engineering,  as  a  part  of  quality  control  exercise,  one  would  like  to 
estimate  the  lifetime  distribution  of  a  day’s  production  of  lightbulbs  from  a 
particular  process.  In  medical  research,  one  of  the  most  common  problems  is  to 
compare  the  lifetime  distributions  of  patients  who  have  undergone  two  different 
treatment  regimens.  One  of  the  regimens  is,  usually,  a  conventional  one  and  the 
other  one  a  new  treatment.  As  an  example,  consider  the  following  problem  we 
have  come  across  recently.  The  conventional  treatment  for  breast  cancer  patients 
has  been  and  continues  to  be  the  one  involving  chemotherapy.  Chemotherapy 
not  only  destroys  cancer  cells,  but  also  good  cells.  The  National  Cancer  Institute 
started  research  into  a  new  treatment:  remove  a  substantial  part  of  bone  marrow 
from  the  patient,  store  it,  and  then  apply  high  doses  of  chemotherapy  until  all 
cancerous  cells  are  destroyed.  The  bone  marrow  is  inserted  back  into  the  bones 
of  the  patient.  The  patient  is  kept  in  an  isolated  ward  until  the  bone  marrow  grows 
back  into  a  sufficient  quantity.  One  of  the  concerns  of  researchers  is  to  compare 
the  lifetime  distribution  of  patients  who  have  undergone  the  new  treatment  with 
that  of  patients  who  have  undergone  the  conventional  treatment. 

Ideally,  in  the  context  of  the  lightbulb  example,  one  takes  a  random  sample 
of  n  bulbs,  sets  them  to  work,  and  follows  them  until  all  of  them  die.  The  observed 

lifetimes  are  usually  written  as  Ti,  T2 .  Tn.  The  next  task  is  to  identify  a 

parametric  distribution  that  fits  the  data  reasonably  well.  A  variety  of  probability 
graph  papers  such  as  Exponential,  Weibull,  Lognormal,  or  Logistic  distributions 
can  be  used  to  plot  the  empirical  distribution  function  of  the  data.  The  empirical 
distribution  function  of  the  data  is  defined  by 

F(t)  =  #  {1  <  i  <  n;  Ti  <  t}/n,  0  <  t  <  «. 

If  the  points  (T,  F(Ti)),  i=1,  2,  ...,  n  lie  more  or  less  on  a  line  in  the  graph,  the 
corresponding  parametric  probability  model  is  fitted  to  the  data. 
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The  full  sampling  plan  of  obsen/ing  each  and  every  unit  in  the  sample  until 
it  dies  is  not  practical  in  many  cases.  Instead,  one  can  select  m  time  points  0  <  ti 
<  ta  <  ...  <  tm  and  inspect  the  units  only  at  each  of  these  m  time  points.  In  such  a 
case,  the  data  consist  of  the  following:  Xi  =  #  units  that  failed  during  (0,  ti], 
Xi  =  #  units  that  failed  during  (tj-i,  ti],  i=2,  3,  ...,  m,  and  Xm+i  =  #  units  that  failed 
during  (tm,  °°).  There  are  some  problems  to  be  resolved  with  this  periodic 
inspection  plan.  (1)  How  many  inspections  m  are  to  be  made?  (2).  How  are  the 
inspection  times  to  be  selected? 

The  first  question  is  addressed  in  Section  2.  The  maximum  likelihood 
estimation  in  multinomial  models  is  discussed  in  Section  3.  In  the  case  when  all 
the  inspection  intervals  are  of  the  same  length  to,  we  will  discuss  an  optimal 
choice  of  to  for  a  given  value  of  m  for  some  parametric  models  of  lifetime 
distribution  in  Section  4.  The  loss  of  information  that  results  in  adopting  the 
periodic  inspection  plan  in  lieu  of  the  full  sampling  plan  is  discussed  in  Section  5. 
The  maximum  likelihood  estimation  for  inten/al-censored  data  is  discussed  in 
Section  6.  In  Section  7,  we  establish  some  results  in  matrix  algebra  that  are 
needed  to  develop  the  linear  model  approach.  Graphical  analysis  for  the  periodic 
inspection  plan  is  discussed  in  Section  8.  In  Section  9,  a  linear  model  approach  is 
developed  to  estimate  the  lifetime  distribution  for  general  periodic  inspection 
plans.  Some  simulation  studies  are  conducted  for  comparing  the  methods  of 
maximum  likelihood  and  linear  model  in  small  samples  in  Section  10.  Finally,  in 
Section  1 1 ,  some  useful  pointers  from  this  study  are  presented. 

Kulldorf  (1961)  is  one  of  the  earliest  researchers  devoted  analysis  of 
interval-censored  data  and  grouped  data.  Kulldorf  (1961)  showed  theexistence 
and  uniqueness  of  maximum  likelihood  estimates  for  interval-censored  data  when 
the  underlying  distribution  is  exponential.  He  also  discussed  the  optimal  length  of 
equi-spaced  inspection  intervals  based  on  the  asymptotic  variance  of  the 
maximum  likelihood  estimator.  Nelson  (1977)  extended  the  work  of  kulldorf 
(1961)  on  exponential  distribution.  See  also  Nelson  (1982)  and  Lawless  (1982). 
Meeker  (1986)  disscused  optimal  length  of  equi-spaced  intervals  in  log  time  for 
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estimating  a  particular  quantile  of  a  Weibull  distribution.  Ostrouchov  and  Meeker 
(1986)  discussed  finite  sample  properties  of  maximum  likelihood  estimators  of 
parameters  of  a  Weibull  distributionin  environment  of  interval-censored  data. 

Kraft  (1992)  discussed  estimation  of  parameters  of  an  exponential 
distribution  when  inspection  times  are  generated  by  a  poisson  process.  Interval- 
cesored  data  appear  in  ecological  and  biological  field  studies.  Some  discussion  of 
methods  is  available  in  journals  devoted  to  this  area.  See  Derleth  and  Sepik 
(1990),  Johnson  (1979),  and  Johnson  and  Christensen  (1986).  For  general 
survays,  see  Sampford  (1952)  and  Deddens  and  Koch  (1988).  For  analysis  of 
interval-censored  data  from  a  nonparametric  angle,  see  Turnbull  (1976),  Pierce, 
Stewart  and  Kopecky  (1979),  and  Akritas  (1988).  In  the  frame  work  proportional 
hazards  and  interval-censored  data  see  Satten  (1996). 


2.  HOW  MANY  INSPECTIONS? 


If  one  wants  to  terminate  the  study  after  a  certain  number  m  of 
inspections,  one  would  like  to  determine  the  minimum  number  of  inspections 
required  in  order  to  estimate  the  parameters  of  the  underlying  distribution.  The 
answer  to  this  question  is  that  the  number  inspections  must  be  at  least  equal  to 
the  number  of  parameters  to  be  estimated.  This  statement  is  amplified  in  the 
following  examples. 

Example  1:  Let  T  be  the  underlying  lifetime  random  variable  with  Logistic 
distribution.  Its  probability  density  function  is  given  by 


fT(t)  = 


ayt^'^ 
(1+  at^  )^ 
0 


t>0, 

otherwise. 


where  a  >  0  and  y  >  0  are  the  unknown  parameters  of  the  model.  Suppose  only 
one  inspection  is  made,  i.e.,  m=1.  The  observable  random  variable  is  then  Xi  =  # 
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units  that  failed  during  (0,  ti].  Note  that  Xi  has  a  Binomial  distribution  with 

(£ 

n  =  sample  size  and  p= — The  likelihood  of  the  data  Xi  =  x  is  given  by 

1+o4| 


L(a,7)  = 

The  likelihood  equations 

and 

lead  to  the  same  equation 


^n^ 


1- 


xr 


1+at|J  [l  +  atj 


“in-x 


— ^nL(a,Y)  =0 
da 


—  £nL(a,Y)  =0 
dy 


dlj 


n  1  +  dtf 

There  is  no  unique  solution. 

Example  2:  Let  T  be  the  underlying  lifetime  random  variable  with  Weibull 
distribution.  Its  probability  density  function  is  given  by 

t>0, 

0  othenwise, 


fT(t)  = 


where  X  >  0  and  y  >  0  are  the  unknown  parameters  of  the  model.  Suppose  only 
one  inspection  is  made,  i.e.,  m=1.  Then  Xi  has  a  Binomial  distribution  with 

-XtY 

n  =  sample  size  and  p  =  1  -  e  K  The  likelihood  of  the  data  Xi  =  x  is  given  by 

iM  = 

The  likelihood  equations 

and 


- , 

1 

a> 

1 

_ 1 

X 

1 

CD 

1 

_ 1 

L  J 

J 

^^nL(X,Y)  =0 
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lead  to  the  same  equation 

There  is  no  unique  solution. 

3.  MAXIMUM  LIKELIHOOD  ESTIMATION  IN  MULTINOMIAL  MODELS 

Every  periodic  inspection  plan  gives  rise  to  multinomial  data.  The 
multinomial  probabilities  are  usually  functions  of  a  finite  number  of  parameters.  It 
is  important  to  identify  the  asymptotic  distribution  of  the  maximum  likelihood 
estimators  of  these  parameters  for  our  work.  The  asymptotic  theory  for  the 
multinomial  distribution  is  somewhat  different  from  the  asymptotic  theory  of  iid 
random  variables.  The  following  result  is  the  bulwark  of  maximum  likelihood 
estimation  in  multinomial  distribution.  The  conditions  of  the  theorem  have 
undergone  some  changes  over  the  years,  but  the  conclusions  are  essentially  the 
same.  See  Fisher  (1928),  Rao  (1958,  1973),  Birch  (1964),  and  Sen  and  Singer 
(1993).  We  report  the  version  given  by  Sen  and  Singer  (1993.  p.  253). 

Suppose  Xi,  Xa,  ....  Xm,  Xm+i  have  a  multinomial  distribution  with 

probabilities  pi(0),  P2(0) . Pm(e),  Pm+i(e),  where 

Xi  +  Xa  +  ...  +  Xm  +  Xm+1  =  n, 

and  each  Pi(0)  is  a  function  of  s  parameters  0i,  0a,  ....  0s,  with  s  <  m+1.  The 
parameter  vector  0=  (01,  0a,  ...,  0s)  of  the  distribution  is  assumed  to  belong  to  a 
known  parameter  space  0,  an  open  subset  of  R®.  Assume  that  each  pj(0)  is 
positive  for  all  0  and  admits  continuous  partial  derivatives  up  to  the  order  2  with 

A 

respect  to  the  components  of  0.  Let  0  be  the  maximum  likelihood  estimator,  i.e., 
0  =  0(Xi,  X2,  ...,  Xm,  Xm+i)  maximizes 


^^nL(?i,Y)  =0 
dy 


n 


6 


nl  -  1 

'-<®>=sTi^n[Pi(9)] 
nxi! 


^.iXi 


i=1 


Let 


""1/2  f) 

aij(e)  =  [Pi(0)]  ^  +  j  =  ‘'-  2.  S. 


Let  A=  A(0)=  (ajj(0)).  Assume  that  A  is  of  full  rank  s.  Let  J(0)  =  nA'A. 


Theorem  3.1  The  asymptotic  distribution  of  (0  -  0)  under  0  e  0  is  the  multivariate 

Normal  distribution  with  mean  vector  6  and  dispersion  matrix  -(A'A)~V 

n 

Some  discussion  is  in  order  on  Theorem  3.1.  Imitating  the  computation  of 
asymptotic  variance  of  the  maximum  likelihood  estimator  in  the  iid  case,  one  can 
compute 

J(0)  =  Eg(-H(0)). 

where 

We  have  observed  that  J(0)  =  H(0).  This  observation  gives  an  alternative  way  of 
obtaining  the  asymptotic  variance-covariance  matrix  of  the  maximum  likelihood 
estimator  for  the  multinomial  case. 


4.  OPTIMAL  CHOICE  OF  THE  LENGTH  OF  EQUI-SPACED 
INSPECTION  INTERVALS 

In  this  section,  we  assume  that  the  inspection  intervals  (0,  ti],  (ti,  h] . 

(tm-1)  tm]  are  of  equal  length  with  common  length  equal  to  to.  We  determine  the 
optimal  value  of  to  based  on  the  asymptotic  variance  of  the  maximum  likelihood 
estimators  of  the  parameters  of  the  model.  Assume  that  the  underlying  lifetime 
distribution  is  Exponential  (0),  0  >  0.  The  probability  density  function  is  given  by 

fx(t;0)  =  0e"^,  t>0. 
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The  likelihood  of  the  data  Xi,  X2, ....  Xm,  Xm+i  is  given  by 


L(0) 


=— fif' 

m+1  HI 

nX:!'=‘' 


nx. 

i=1 


]  [' 


iX 


m+1 


The  maximum  likelihood  estimator  0^  of  0  is  given  by 


*0 


0n  = — ^n| 


m 


mn+  S(i-m)Xi 

_ M _ 

m 

mn+  X(i-1-m)Xj 
V  i=1 


with  asymptotic  variance 


Var(0n) 


(e^^o  - 1)^ 
nt2e®^o(1-e'"^®*o)' 


The  formulas  presented  here  match  with  those  of  Nelson  (1977).  For  a  given 
value  of  6  and  m,  nVar(0n)  is  minimized  with  respect  to  to.  nelson  (1977) 
provides  optimal  values  of  0to.  We  tabulate  the  optimal  values  of  to  for  set  of  0 
and  m  values  (Table  4.1).  The  table  also  lists  nVar(0n)  under  the  column  VAR. 
Our  table  is  more  useful  in  answering  several  questions  from  a  practical 
perspective  than  Nelson  table.  If  one  decides  how  many  visits  m  to  be  made  and 
has  arough  idea  about  the  value  of  0,  the  length  to  of  the  inspection  intervalscan 
be  read  off  the  table.  If  one  knows  the  value  of  0  to  be  within  20%  of  some 
specified  value  0o,  one  could  use  the  optimal  length  to  associated  with  0o  in  this 
sampling  plan  without  increasing  the  asymptotic  variance  substantially.  We 
illustrate  this  point  with  a  few  examipes. 
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Table  4.1 

Optimum  Inspection  Interval  Length  to  For  Different  Values  of  0  and  m 
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We  now  take  up  the  case  of  the  Weibull  distribution.  The  probability 
density  function  is  given  by 


fT(t) 


0 


ie-xt^ 


t>0, 

otherwise, 


Likelihood  of  the  data: 


-  n*  Ji- 

n 

nxi!'-' 

i=l 

Loglikeiihood  of  the  data: 


Xi 


(n-SXj) 


i=l 


m 


^(0)  =  const.  +  XXi  ^n| 
i=l 


.-A(i-l)t5  _  -Aitg 


ArntJXm+i. 


First-order  derivatives  of  the  loglikeiihood: 


dm_  ™  itje 


dX 


=  I 


i=l 


Xj  -mtJX^+l, 


and 


sue)  ™  ^ . 

.  -;u.5 - X;-Amt5(<„.o)X..„ 


These  are  two  non-iinear  equations  in  X  and  y  which  can  be  solved  by  using  an 
iterative  procedure  such  as  Newton-Raphson  method.  It  can  be  shown  that 
the  information  matrix  is  given  by 


1(0)  = 


Si2 

$22 


where 
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Si  1  =  -E  — =  nY  J — , 

_nx,2y  ito-[('-1)toF^n  [(i-1)to]]^ 

i=i  LMito)^  _eM('-'')‘o]] 


and 


32«g)  >»[(i'o)’'<''(ito)-[(i-1)torM0-1)<o]T0'o)’'-[(i-1)tof 

S-|  2  =  ”E  -r - =  nX  Y  ‘^= - P - =L - - 

=  ®ai- 

The  inverse  of  the  information  matrix  is  given  by 


For  given  values  of  m,  A,,  and  y,  the  trace  of  the  information  matrix,  apart  from  n, 
is  minimized  with  respect  to  to.  The  optimum  values  are  tabulated  in  Table  4.2, 
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0.333  0.843 _ 10  0.71  2.4615  0.61  1.9320  0.52  1.7333  0.47  1.6237  0.43  1.5523 

0.500  0.447 _ 10  1.37  2.6757  1.19  1.8026  1.08  1.4917  1.00  1.3262  0.94  1.2226 

1.000  0.100  10  7.74  3.4273  6.53  2.0478  5.68  1.5928  5.06  1.3663  4.56  1.2304 

1.500  0.027 _ 10  8.53  4.5951  6.93  2.8490  5.84  2.3003  5.07  2.0376  4.48  1.8853 

2.000  0.008  I  10  I  8.43  5.9753  I  6.6  4.0450  5.41  3.4257  4.60  3.1268  4.00  2.9539 


5.  LOSS  OF  INFORMATION 

One  of  the  questions  that  arises  in  the  case  of  periodic  inspection  plan  is 
that  how  much  information  we  lose  if  we  adopt  this  plan  in  lieu  of  the  full  sampling 
plan.  We  proceed  to  answer  this  question  in  the  following  way. 

1 .  Evaluate  the  asymptotic  variances  of  the  maximum  likelihood  estimators  of  the 
parameters  of  the  underlying  model  under  each  of  the  full  sampling  and 
periodic  inspection  plans. 

2.  Compare  the  variances  of  the  plans. 

A  comparison  is  made  by  computing  the  efficiency  of  the  full  sampling 
relative  to  the  periodic  inspection  plan.  We  consider  only  the  periodic  inspection 
plan  in  which  the  inspection  intervals  are  equi-spaced.  The  relative  efficiency 
(RE)  is  defined  as  the  ratio  of  the  variance  of  the  estimator  under  the  full  sampling 
plan  to  that  under  the  periodic  inspection  plan  in  the  single  parameter  case  under 
the  optimal  choice  of  to.  It  is  defined  as  the  ratio  of  the  trace  of  the  asymptotic 
variance-covariance  matrix  under  the  full  sampling  plan  to  that  under  the  periodic 
plan  again  under  the  optimal  choice  of  to  in  the  multi-parameter  case. 

First,  we  take  up  the  case  of  the  Exponential  distribution  with  parameter  0  > 
0.  Let  m  stand  for  the  number  inspections  and  to  the  corresponding  optimal 
common  length  of  the  inspection  intervals.  For  each  value  of  m,  we  will  compute 
the  efficiency  of  the  full  sampling  plan  relative  to  the  periodic  inspection  plan. 
Before  this,  we  need  to  determine  the  maximum  likelihood  estimator  of  0  under 
the  full  inspection  sampling  plan.  Under  the  full  sampling  plan,  let  Ti,  Ta, ...,  Tn  be 
the  lifetimes  of  the  units.  The  maximum  likelihood  estimator  0n  of  0  is  given  by 

0  -  " 

0n-  — . 

STi 

i=1 

0^ 

and  the  asymptotic  variance  of  ©p  is  given  by  — . 

Therefore,  the  efficiency  of  the  full  sampling  plan  relative  to  the  periodic  inspection 
plan  is  given  by 
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RE=  °  et  2 - 

(e®»o  _i)2 

Table  5.1  shows  the  optimum  value  of  to  and  the  corresponding  relative  efficiency 
for  e  =  0.1 , 0.2, ...,  0.9, 1 , 2,  ...,10  and  for  m  =  1 , 2,  3,  4,  5,  6. 
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Table  5.1 
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From  Table  5.1,  we  notice  that  the  relative  efficiency  increases  as  m,  the 
number  of  inspections,  increases.  Also,  the  relative  efficiency  increased  by  1 6%  if 
two  inspections  are  made  rather  than  one  inspection.  This  increment  becomes 
smaller  as  m  increases.  We  notice  also  that  optimum  value  of  to  is  approximately 
equal  to  the  expected  value  of  the  lifetime  distribution  T  when  three  inspections 
are  made.  So,  if  the  lifetime  distribution  is  Exponential,  one  can  make  three 
inspections  at  times  to,  2to,  and  3to,  where  to  is  approximately  the  expected  value 
of  the  lifetime  distribution. 

We  now  take  up  the  case  of  the  Weibull  distribution.  Consider  an 
inspection  plan  consisting  of  m  inspection  times  ti  <  ta  <  ...  <  U.  Let 

X|  =  #  units  that  failed  during  the  time  interval  [ti.i,  tj),  i=1 , 2, ...,  m, 
and 

Xm+1  =  #  units  that  failed  during  the  time  interval  (tm,  “o), 
with  the  convention  that  (tn,  tj]  =  (0,  ti],  for  i  =  1.  The  parameter  vector  is 
0  =  (X,,y),  and  the  parameter  space  0  is  identified  as 

R+  X  =  je  =  (X, y):  X  >  0,  y  >  o}.  it  is  dear  that  Xi,  Xa,  ...,  Xm,  Xm+i  have  a 
multinomial  distribution  with  probabilities 

Pi(e)  =  e"^*"i-e"^'",  i  =  1,  2 . m, 

and 

Pm*l(S)  =  e-"". 

All  the  conditions  of  Theorem  4.3.1  are  met  in  this  case.  The  likelihood-based 
information  matrix  1(9)  is  computed  after  the  following  derivation,  with  the 
convention  that  0£n0  =  0. 

Likelihood  of  the  data: 
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L(e)  = 


n! 


m 

n 


m+1  . 

nxii"" 

i=1 


e 


nXi 


e 


m 

>-IXi) 
1=1 


Loglikelihood  of  the  data: 

^(0)  =  const.  + 


fx, 


-XilX 


m^m+1- 


First-order  derivatives  of  the  loglikelihood: 


d£{Q) 

dX 


m  tj^e-^*'  -t^  e-^*>-i 
=  I  - ..  ■■ — Xi-t^,X 


=1 


e 


and 


These  are  two  non-linear  equations  in  X  and  y  which  can  be  solved  by  using  an 
iterative  procedure  such  as  Newton-Raphson  method.  To  obtain  the  Fisher-Rao- 
Birch  information  matrix  J(0)  in  this  two-parameter  case,  we  need  to  identify  the 
entries, 


if  j  =  1, 
ifj  =  2, 


1 

apj(0) 

VpiW 

< 

di 

1 

3pi(g) 

.VPiW 

dy 

i  =  1, 2 . m,  m+1. 

Therefore,  the  matrix  A  is  given  by 


A  =  A(0)  = 


aii  ai2 

^21  ®22 


.^m+11  ®m+12j 
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where 


A'A  = 


m+1  „  m+1 

I  a2  Xa,ai2 


Xaiiai2  X  a| 


J(§)  =  nA'A  = 


n  X  n  Xajiaig 
i=1  i=1 

m+1  m+1 

n  Zaiia|2  n  X  afp 


Jll(^)  Jl2(^) 
J2i(0)  J22(®) 


Jll(e)  -  n  X  aj^  -  X;^ 

i=1  i=iPi(0)V  j 

‘nf—  ^  .,y»l 

1=1  e  '^*1-1  -e  ^*1  ^ 

=  nX^= - =L+-ntm^e“^*^, 


J.a(5)  =  ns’a|  =  I^^f 

i=1  i=lPi(0)v  ^  J 

m+l  H  r  ^.y  *  Y 

=  >^1  — }^ - TTT  -^tM(^n^-i)®  +A,tT(Mj)e 

i=i  e  ''■^1-1  ^ 


=  nx25; 
i=l 


t7(£nti)e"^*'^  -  tj’[.^(Mj_i)e"^^'^-i 


•  +  nA,2tS^(Mm)e 
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m+1  m 


J,2(e)  =  n2ai,ai2=X-U^ 

i=1  i=lPi(0)l^  oX 


f  t7\\\ 


A  ^  ; 


m+1 

=  nl 

i=1 


e-wJ-i  -e-"/ 


m 

=  nXIi 

i=1 


ire"*-'!'  tr(«nti)e‘"i'  -t? 


A  special  case  of  this  result  is  of  interest.  Suppose  all  inspection  intervals 
have  the  same  length  to,  i.e.,  ti  =  ito,  i  =1,  2,  ....  m.  It  can  be  shown  that  the 
information  matrix  is  given  by 

Sii  Si2 

.S21  S22. 


1(0)  = 


where 


1  =  -c  — 5—  =  n  X  - - 


y  * 


e  '''v-u/ 

_  r 

S22  =  “E - ;r“  =  X  - - f= - =5 - 

^  i=1  _qM('“^)*o] 


and 


axdy 


i=1 


(ito)^^n  (ito)-[(i 

i  -  1)lof  £n[(i  -  1)to]J(ito)''  -  [(i  -  1)tof  ’ 

e^ito)^  _e^[(*-'')‘or 

=  $21. 

The  inverse  of  the  information  matrix  is  given  by 
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Sii 

Sl2' 

S21 

S22. 

We  now  obtain  the  information  matrix  of  the  maximum  likelihood  estimators 
of  the  parameters  of  the  Weibull  distribution  under  the  full  sampling  plan.  Let  the 
lifetime  T  have  a  Weibull  distribution.  The  probability  density  of  T  is  given  by 


fT(t)  = 


0 


t>0, 

otherwise, 


where  A.  >  0  and  y  >  0  are  the  unknown  parameters  of  the  model. 

Under  the  full  sampling  plan,  select  n  units  at  random.  Switch  them  on. 

Watch  them  until  all  of  them  fail.  Let  Ti,  Tz . Tn  be  the  lifetimes  of  the  units  1 , 2, 

...,  n  and  h,  ta, ....  tn  the  observed  failure  times.  The  likelihood  function  is  given  by 


L{X,y)  =  {Xyf  Hti 
J=1  . 

^(X,Y)  =  ^nL(X,Y) 


e  ‘=1  . 


n  n 

=  n^nA,  +  n^ny  +  (y  -  1)X  • 

i=1  i=1 


Differentiate  the  loglikelihood  with  respect  to  X  and  y,  and  equate  the  derivatives  to 
zero: 


a^(X,y) 

ax 

a^(x,y) 

dy 


2t?'=o, 


i=1 


n  n  n 

=  -  +  X^ntj-XXtj'^Mj  =0. 
y  i=1  i=1 


These  are  two  non-linear  equations  in  X  and  y  which  can  be  solved  by  using  an 
iterative  procedure,  such  as  Newton-Raphson  algorithm. 

The  second-order  partial  derivatives  of  the  loglikelihood  are  given  by 


a^^(X,y) 

ax^ 

a^^(X,y) 

ay^ 


--xhnientf, 

Y'"  i=1 
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and 


d^t{X.y) 

3»r 


The  negative  of  the  expectation  of  the  second-order  partial  derivatives  are  needed 
to  be  evaluated.  The  simplest  one  is  given  by 

d'^£iX,y)  n 

■ 


The  other  expectations  require  additional  work.  First,  we  need  to  introduce  the 
Digamma  function  (Abramowitz  and  Stegun,  1970). 

For  any  z  >  0,  the  Gamma  function  r(z)  is  defined  by 


the  Digamma  function  by 


r(z)  =  J  t^' V*dt, 

0 


'P(z)  =  —  ^nr(z) 
dz 

=  -^^n?t^"V‘dt 
dz  i 


J(M)t^“''e"*dt 


_  0 


r(z) 


and  the  derivative  of  the  Digamma  function  by 

T'Cz)  =  ^'F(z)  =  J(M)  e~*dt.^^^^‘^^  ^^nt)-t^  V^z) 
0  (r(z))2 


=  J(M) 


2  t^-^e- 


■<lt-H'(z)j^-^(M)dt 


r(z) 


r(z) 


=  J(M)2ll^t-'F2(z). 

n  ^ 


We  now  compute 
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-g(3^<(X,7))  ^  n  ^ 1  ^  n  ^ J (<n,)2,Y j_^T-ie-M''-’^ 


ay^  y-  y 

Let  x  =  Xt'^.  Then  dx=  ^Tt'^^'^dt,  and 


0 


-i2 


E(Tl'(«nT)2)  =  J 

oL  ^ 


X  -X  j 

—  e  *dx 


oo  . 

=  J— (£n  X  -  ^n  X)^  j  e'^dx 


OY 

_l_ 

Xy‘ 


X 

j  (^n  x)^x  e'^dx-2^n  A,j(£n  x)  x  e’^dx  +  (£n  ?i)^  Jxe'^dx 


10 


=  ^|'P'(2)  +  ^^{2)-2^{2)inX  +  (^n  X)2|. 


Therefore, 


2 

+  vp2(2)-2'F(2)^nX  +  (^n>.)2|. 


Further, 


dy2 


2 

=  nE(TY(£nT))  =  nj  tY(^n  t)XYtT^-‘'e-^*"dt. 
dXdy  J 

r-1. 


Let  x  =  XtY.  Then,  dx=  ^TtY'^dt,  and 


n  ^ 

=  —  J  (^n  X  -  ^n  X)  xe'^^dx 
^Y  0 


n 

~Xy^ 
n 


J  (£n  x)  xe’^dx  -  J  (^n  X)  xe’^'dx 
0  0 


=  :^{'F(2)-^n;.}, 

Xy 

Therefore,  the  information  matrix  is  given  by 
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I  =  n 


1 


2 

Xy 


('F(2)  -  £nX) 


^(4^(2)  -  £n;i)  +  ^'(2)  +  (4^(2)  -  inXf'^ 


Xy  Y* 

The  determinant  A  of  the  information  matrix  is  given  by 

A  =  +  '*"(2)  +  ('*'(2)  -  -  ^('•'(2)  -  fr>xf 

^  X,  Y 


X‘^y‘ 
n 


xV'  '  xV  ’ 

The  inverse  of  the  information  matrix  is  given  by 


(1  +  'F'(2)) 


n4"(1) 


n 


(y(2)-^nX)^'|  (4^(2) -^n>.) 


Xy 


'F'O) 
(4^(2) -£nX) 
^'(1) 


4>'(1) 

.2 


^'(1) 


^11  2.12 
L^21  ^22  ^ 


say. 


The  efficiency  of  the  full  sampling  plan  relative  to  the  periodic  inspection  plan  is 
defined  by 

^11+^22 

Sii  +S22 

Table  5.2  shows  the  optimum  value  of  to  and  the  corresponding  relative  efficiency 
for  different  combinations  of  y,  X  and  m. 
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Table  5.2 

Optimum  Inspection  Interval  Length  to  and  Relative  (RE)  of  Full  Samling  Plan 
Versus  Periodic  Inspection  Plan  for  different  combination  of  X,  y  and  m. 


CM 

CD 

N 

1 

1 

00 

CD 

00 

CO 

o> 

in 

CO 

CO 

00 

in 

o 

pi» 

CO 

CD 

oo 

CD 

CO 

CM 

CM 

CD 

00 

CO 

■ 

iy 

CM 

ay 

CO 

CO 

00 

CM 

h- 

SO 

in 

■ 

Q 

00 

in 

in 

00 

II  CE 
c 

o 

o 

o 

O 

CO 

CD 

o 

1 

|| 

d 

d 

d 

d 

d 

d 

w 

CM 

o 

CM 

CO 

CD 

in 

■ 

o 

CO 

CO 

00 

o 

o 

■M* 

CO 

in 

■ 

Ea 

o 

CD 

in 

o 

o 

o 

o 

o 

o 

o 

o 

■ 

cvi 

o 

d 

B 

B 

Q 

CD 

CD 

CO 

CD 

CO 

CO 

o 

CO 

CO 

CO 

CM 

o 

LXJ 

h- 

in 

CO 

O 

00 

T* 

in 

00 

CO 

00 

CO 

CO 

in 

in 

00 

CVJ 

ay 

CO 

in 

CM 

CD 

CD 

h- 

f: 

o 

00 

f:: 

m  cr 

h- 

CO 

CO 

00 

CO 

■M- 

CO 

m 

CO 

II 

c 

O 

o 

o 

o 

o 

O 

O 

o 

d 

d 

d 

d 

d 

d 

o 

in 

CD 

in 

in 

in 

Tf 

CO 

h- 

o 

CO 

o 

■M* 

in 

'M’ 

in 

CO 

Tj* 

o 

o 

o 

CO 

o 

o 

o 

o 

o 

o 

o 

CM 

CM 

CM 

o 

in 

B 

s 

N 

CO 

CO 

00 

00 

N 

CO 

CO 

T“ 

00 

CO 

CM 

LU 

N 

in 

tT 

CO 

CO 

CD 

CD 

CD 

O 

h- 

00 

CO 

00 

O 

ay 

in 

CM 

CO 

'M- 

00 

Tf 

CM 

CD 

CO 

in 

T— 

CD 

CO 

CO 

CO 

CO 

r- 

in 

Tj- 

in 

h- 

in 

CO 

Tf 

in 

h- 

n 

c 

o 

6 

o 

o 

o 

o 

o 

O 

o 

o 

O 

d 

o 

o 

d 

o 

o 

CJ) 

CO 

CO 

in 

CO 

CM 

T“ 

CM 

00 

00 

in 

in 

in 

in 

in 

CO 

N 

CD 

in 

o 

CO 

00 

o 

o 

o 

o 

o 

o 

o 

CM 

CM 

cvi 

o 

in 

Q 

Q 

CO 

in 

CO 

in 

CM 

CO 

CD 

CD 

00 

CD 

LU 

CM 

CM 

<yy 

O) 

CO 

CD 

o 

CO 

CO 

CO 

CM 

m 

CD 

CM 

h- 

CO 

00 

CO 

oo 

o 

00 

CD 

CM 

00 

o 

CO  CO 

CO 

in 

in 

in 

CO 

in 

CO 

CO 

CO 

CM 

CO 

CO 

II 

p 

6 

CD 

CD 

CD 

CD 

CD 

d 

o 

O 

d 

d 

d 

o 

O 

d 

CO 

CO 

in 

00 

in 

CM 

CO 

T- 

CD 

CO 

CO 

in 

in 

CO 

CO 

CD 

in 

CO 

T- 

CO 

co 

T— 

in 

CD 

CO 

CO 

CD 

CD 

CD 

CD 

CD 

d 

CO 

d 

CO 

CO 

cb 

CM 

ay 

CO 

CD 

h- 

CO 

N. 

CO 

CM 

o 

in 

Tf 

CM 

LU 

TJ- 

ay 

CM 

00 

CD 

CO 

O 

00 

o 

CO 

CD 

h- 

CO 

CO 

CD 

in 

CO 

o 

O 

o 

o 

ay 

CD 

O 

CVJ  CC 

in 

Tt 

CO 

CO 

Tf 

rf 

CM 

CO 

•M- 

CO 

CM 

CM 

II 

d 

CD 

CD 

CD 

CD 

d 

o 

O 

d 

d 

d 

d 

d 

d 

& 

£ 

cn 

CO 

CO 

CO 

in 

o 

in 

T- 

CO 

CO 

CO 

CO 

1^ 

00 

00 

h- 

CM 

IcM 

h- 

q 

r^ 

in 

CD 

CD 

d 

CD 

CD 

d 

d 

m 

B 

B 

d 

N 

00 

cb 

p 

uT 

in 

in 

in 

in 

in 

10 

o 

o 

o 

T“ 

o 

m 

o 

CD 

in 

CO 

CM 

o 

f: 

CO 

N 

o 

00 

ir“ 

o 

in 

00 

CO 

CO 

o 

CO 

Tf 

o 

CM 

CD 

CO 

■M- 

o 

00 

q 

CO 

CM 

o 

o 

00 

O 

o 

■ 

T- 

T~ 

o 

o 

T“ 

d 

o 

o 

o 

o 

d 

d 

O 

o 

CO 

o 

o 

o 

o 

CO 

o 

o 

o 

o 

CO 

o 

o 

o 

o 

1 1 

CO 

o 

o 

o 

o 

CO 

o 

o 

o 

o 

CO 

o 

o 

o 

o 

CO 

in 

q 

q 

o 

CO 

in 

q 

in 

o 

CO 

in 

q 

q 

o 

mi 

CD 

CD 

T“ 

cvi 

d 

d 

CM 

o 

d 

1" 

1F“ 

cvi 

CN 


Comments: 


Relative  efficiency  increases  with  increasing  values  of  m,  the  number  of 
inspections.  This  is  to  be  expected.  The  relative  efficiency  does  not  exceed  0.86 
in  ail  cases  considered.  If  Var(T)  is  large,  relative  efficiency  is  poor.  The  relative 
efficiency  is  very  much  stable  for  m  =  4,  5,  and  6.  Four  inspections  seems  to  be 
ideal.  The  optimum  length  to  decreases  with  increasing  values  of  m.  This 
phenomenon  is  natural  to  expect. 


6.  MAXIMUM  LIKELIHOOD  ESTIMATION  FOR  INTERVAL-CENSORED  DATA 

In  this  section,  we  discuss  the  method  of  maximum  likelihood  that  is 
normally  used  to  estimate  the  lifetime  distribution  parameters  for  interval-censored 
data.  In  most  cases,  the  likelihood  equations  do  not  have  explicit  solutions.  We 
record  the  likelihood  equation  and  asymptotic  variances  of  estimators  for  each  of 
the  six  parametric  distributions:  Exponential,  Weibull,  Lognormal,  Logistic,  pareto, 
and  the  Extreme  value.  If  the  sample  size  is  reasonably  large,  one  could  use  the 
method  of  maximum  likelihood.  The  equations  are  and  formulas  are  in  place. 

In  Section  8,  we  introduce  an  alternative  method,  the  Linear  Model 
Approach  for  inten/al-censored  data,  in  which  an  explicit  formula  is  obtainable  in 
the  estimation  of  the  lifetime  distribution  parameters.  We  will  show  that  the  new 
method  give  consistent  estimators.  Recall  that  the  inspection  plan  consists  of  m 
inspection  times  ti  <  ta  <  ...  <  U,  Xi  =  #  units  that  failed  during  the  time  interval 
(ti.i,  tj],  i=1, 2 . m,  and  Xm+i  =  #  units  that  failed  during  the  time  interval  (tm,  «>). 


6.1  Maximum  Likelihood  Estimation  for  the  Exponential  Distribution 

The  maximum  likelihood  estimator  Sp  of  0,  where  0  is  the  parameter  of  the 
underlying  Exponential  distribution,  is  a  solution  of  the  likelihood  equation. 


-ti_ie‘®*i-i  +tie"®* 


j 


m 


d0 


i=1 


-  (n-lXi)tp,=0.  (6.1.1) 

i=1 
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The  equation  is  non-linear,  which  can  be  solved  using  an  iterative  procedure.  The 
asymptotic  variance  of  0^,  Using  Theorem  3.1,  is  given  by 


ASYVAR(en) 


1 


1 


-i-nt^e 


(6.1.2) 


In  the  special  case  where  tj  =  ito,  (6.1.1)  has  an  explicit  solution  given  by 


-  1 
0  =  — ^n| 

^0 


m 

mn+  X(i-m)Xi 

i=1 _ 

m 

mnn-  Ki-1-m)Xj 
i=1  J 


and  the  asymptotic  variance  of  0n  in  this  case  is  given  by 


ASYVAR(0n)  = 


(e^^'o  _i): 


nt2e®*o(1-e‘'^®*o)' 


6.2  Maximum  Likelihood  Estimation  for  the  Weibull  Distribution 

The  maximum  likelihood  estimators  X  and  y  of  X,  and  y,  respectively,  where 

X  and  y  are  the  parameters  of  the  underlying  Weibull  distribution,  are  a  solution  of 
the  likelihood  equations. 


9^(0) 

ax. 


tfe 


e 


-M?-, 


^m^m+1  “  0, 


and 


^  i=i  i=i  - 1 


0. 


These  are  two  non-linear  equations  in  X  and  y  which  can  be  solved  by  using  an 
iterative  procedure  such  as  Newton-Raphson  method.  Even  in  the  special  case  ti 
=  ito,  there  is  no  explicit  solution.  The  asymptotic  variance-covariance  matrix  is 
given  by  the  inverse  of  the  matrix,  using  Theorem  3.1 , 
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nA'A  = 


Sii  Si2 

1^21  S22 


where 


m 

S11  =  nX 

i=i 


i2 


e 


m 


_ d _ L  nt  6 

1  ^  ''‘m  ®  > 


$22  =  nA,  ^ 
i=1 


t 


r 

(<nt|)-tr,{^nt|.i)]‘ 


e"?  -e^’- 


and 


i=1 


e"? 


6.3  Maximum  Likelihood  Estimation  for  the  Lognormal  Distribution 

The  probability  density  function  of  Lognormal  distribution  is  given  by 


fT(t)  = 


1 


taV^ 

0 


e 


2o 


t>0, 

otherwise. 


The  parameter  vector  is  0  =  (|x,a^)  and  the  parameter  space  ©  is  identified  as 
=|0  =  ()a,a^):  -oo<jx<oo,  a>o|.  The  random  variables  Xi,  X2,  Xm, 
Xm+1  have  a  multinomial  distribution  with  probabilities, 

Pj(e)  =  <D  — _<i>  — LJ_Hi  I  1  =  1,  2. 

V  j  V  ) 

and 


m, 


where 
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Likelihood  of  the  data: 


<I>(x)  =  j 


^  e  2  dz. 


~  n*  JIl 

'-(9)=di-n 

nxi!'=' 

1=1 


Loglikelihood  of  the  data: 
^(0)  =  const.  + 

i=1 


-<D 

■Xi 

1-<D 

=1 

1 

E 

c 

.  V  CT  J 

1  ^  ). 

m 

(n-IXi) 

i=i 


1-0) 


V 


First-order  derivatives  of  the  loglikelihood: 


d£(G)  ^  g 


(l>(^ntm) 


i=i  o 


Un\.\  -p"' 

_<D 

r^nti_i-p^ 

-/Vj  T 

1-0) 

^^ntm 

<7  J 

1  <7  J 

1  cr  J 

^  )) 
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30=“  ^ 


2a‘ 


(Mm  -4)(})(Mm) 
2<y^ 
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=1 

1 

C 

-0) 

/Nj  T 

1-0) 

1 

E 

c 

V  J 

I  J 

i  (J  ; 

X, 


_ 

where  m)  =  — ^e  2a^ 
av27t 


(t-n)' 


The  maximum  likelihood  estimators  ji  and  a^of  |i  and  a^,  respectively,  are  found 

by  equating  the  first-order  derivatives  to  zero.  The  resulting  equations  are  non¬ 
linear  which  can  be  solved  by  using  an  iterative  procedure  such  as  Newton- 
Raphson  method.  The  asymptotic  variance-covariance  matrix  is  given  by  the 
inverse  of  the  matrix, 


nA'A  = 


Sii  Si2 

$21  S22 


where 
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m+1  -i 

S11  =  I  ^ 
i=i  Pi(0) 


0Pi(e)^^ 
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i=i  Pi(0) 


and 
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1-0 

:  1 

E 

c 

1  J 
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I 
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1  20^^  20*^ 
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Mj  -p"' 
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-0 

1  J 

Mrn,-:P)^2(^nt,„) 
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1-0 


6.4  Maximum  Likelihood  Estimation  for  the  Logistic  Distribution 
The  probability  density  function  of  Logistic  distribution  is  given  by 


fT(t)  = 


t>0, 


(1  +  A.tT^)2 
0  otherwise. 
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The  parameter  vector  is  6  =  (A,, 7),  and  the  parameter  space  0  is  identified  as 
R+ xR'^  =|0  =  (A„y):  ^>0,  a>o}.  The  random  variables  Xi,  X2,  ....  Xm,  X, 
have  a  multinomial  distribution  with  probabilities, 


m+1 


Pi(0)  =  P(tj_i  <T<tj), 
1  1 


^+x\y  ^+uy. 

I  I- 1 


.  i  =  1,  2 . m, 


and 


Pm+l(®)  ~  ^m) 

1 


Likelihood  of  the  data: 


ni 

nxi!'=’ 

i=1 


— fi 

m+1  A  A 

nxii'=' 

i=1 


1  1 

Xi 

( 

1 

1  +  A,tj^^  1+X.tj^ 

Xi 

xxj  -  Xtl^ 

1 

(i+uj)(unl,)_ 

m 

-S) 

i=1 


_^m+1 


Loglikelihood  of  the  data: 

^(e)  =  const.  +  “X  +  +  + 

i=l  '■  ■'  i=i  '  J 


-Xm+i^n(1  +  ^n(XtJ,)). 
First-order  derivatives  of  the  loglikelihood: 


de{Q)  ^  ^ 
^=1 


1  t/ 


t7 

‘i-1 


^  ^+x\y  ^+uy_^ 


t^ 

Xj  — ini_x 


m+1' 


and 
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3^(0)  m 

ay  “it,  .^.y,) 

A,tm(^ntm) 


i=1 


1  +  X\l 


X 


m+V 


I  +  Xtj  1  +  Xt7.^ 


The  maximum  likelihood  estimators  X  and  y  of  X  and  y,  respectively,  are  found 

by  equating  the  first-order  derivatives  to  zero.  The  resulting  equations  are  non¬ 
linear  which  can  be  solved  by  using  an  iterative  procedure  such  as  Newton- 
Raphson  method.  The  asymptotic  variance-covariance  matrix  is  given  by  the 
inverse  of  the  matrix, 


nA'A  = 


Si2 

S22; 


where 


m+1 


1 


S  Pi  (5) 


aPi(9)''^ 

dX 


m 

=  I 

i=1 


i2 


t^ 

N-1 


*?' 


i+xtj  i+x\y_^ 


■  + 


^2y 


1  +  ^tm 


3  ’ 
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m+1 


1 


S22  =  X  ,3. 

i=1  Pi(0) 


^Pi(9) 

ay 


^.n2 


m 

=  I 

i=1 


l2 


uxtl, 


1+ 


xxj] 


1 


1 


1  +  xtj^ 


j.  ^  ^rn 

vi3  ’ 
1  +  XtJ, 


and 


o  _o  1  api(0)  api(0) 

^12  -^21  =  L  .  — 

i=i  Pi(0)  dy 


m 

=  I 

i=1 


{y  t^ 

‘i-1  N 

XXjinXi 

1  +  A.tj^]  ^ 

[uuyf 

1  1 

6.5  Maximum  Likelihood  Estimation  for  the  Pareto  Distribution 

The  probability  density  function  of  Pareto  distribution  is  given  by 


frO)  = 


ya^ 

'  tY+1 
0 


t>a, 

otherwise, 


where  7  >  0  is  the  unknown  parameter  of  the  model.  Assume  that  a  >  0  is  known. 
In  this  case,  the  first  inspection  interval  is  (a,  ti].  The  random  variables  Xi,  X2,  ..., 
Xm,  Xm+1  have  a  multinomial  distribution  with  probabilities, 

Pi=P(ti_i<T<ti), 
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Likelihood  of  the  data: 


Pm+1  “  >  ^m) 


i=1 


Loglikelihood  of  the  data: 

\f  a  V  raVl  r  a  1 

^(Y)  =  const.  + 

j=1  Vm-1/  vm/  L‘m. 

First-order  derivatives  of  the  loglikelihood: 


The  maximum  likelihood  estimator  y  of  y,  is  found  by  equating  the  first-order 

derivative  to  zero.  The  resulting  equation  is  non-linear,  which  can  be  solved  by 
using  an  iterative  procedure  such  as  Newton-Raphson  method.  The  asymptotic 
variance  of  y  is  given  by  the  reciprocal  of  the  information  number  l(y), 


6.6  Maximum  Likelihood  Estimation  for  the  Extreme  Value  Distribution 
The  probability  density  function  of  Extreme  value  distribution  is  given  by 
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The  parameter  vector  is  0  =  (0,y),  and  the  parameter  space  0  is  identified  as 
R  X  =  (0, y):  -  oo  <  0  <  oo,  y  >  oj.  Let  T  =  e^ .  It  can  be  shown  that  the 

probability  density  function  of  T  is  given  by 

where  X  =  e”^,  which  is  the  probability  density  function  of  the  Weibull 
distribution  with  parameters  y  and  X.  Therefore,  to  find  the  maximum  likelihood 
estimators  for  the  parameters  of  the  Extreme  value  distribution,  we  fit  the  Weibull 

distribution  to  the  T=e^  variable.  The  estimator  y  of  y  is  equal  toy, and  the 
estimator  0  of  0  is  equal  to  {-inX)/y.  By  using  the  A-method,  the  asymptotic 
variances  of  y  and  0  are  given  by 

asy  var(y)  =  S22 


and 


asvar(0) 


2..2 


X^^y 


_  1  q/  .  £nX 

--Tr-^Oii+ - ;5— 022  -2 — 


Xy'^ 


12. 


where 


l(0)  = 


■Sii  S^2 

S22 


is  the  asymptotic  variance-covariance  of  X  and  y. 
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7.  SOME  RESULTS  FROM  MATRIX  ALGEBRA 

In  this  section,  we  establish  some  results  in  matrix  algebra.  These  results 
will  be  used  in  the  subsequent  sections.  Let  the  mxm  symmetric  matrix  2  have 
the  following  special  structure: 

0-|  0-|  6-(  0-|  .  .  .  0.|  0.| 

02  02  ^2  ®2 
03  •  •  •  03  03 

04  04  .  .  .  04  04 

2  =  1  05  ...  05  05 


0(m-1)  0(m-1) 

0 


m  ) 


(i)  The  determinant  A  of  this  matrix  is  given  by 


m 


n(0i-0i-i). 

i=1 

with  the  convention  that  0o  =  0. 

ii)  The  inverse  matrix,  =  (ay)  of  the  matrix  Z  is  a  tri-diagonal  matrix  and  is 
given  by 

(0i+1  ~  0i-i) 


an  = 


(0i-0i_i)(0M-0i) 


,  i  =  1,  2,  ...,  m-1. 


^mm  “ 


(0m  -0m-i) 


and 


^i,i+1  -  - 


1 


(0i+1-0i) 


Let 
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‘tr 

’yi' 

ti 

yi 

t.  = 

• 

and  y  = 

. 

-^m_ 

.Vm. 

be  mx1  vectors.  Then 


m 

I 

i=1 


m 


■*'2X 

i=2 


and 


m  m 

(t*Z  Y)  =  ^  tjyjOjj  +  (tj-iyj  +  tjyj_i)aj_-|  j. 
i=l  1=2 

(iv)  Let 

‘1  U' 

1  t2 


1  t 


m 


Then 


(t'2:"‘‘t)  = 


m  m  m  m 

X  Cii  +2Xaj_i,i  X  +  X  (ti-1 +ti)<7i-i,i) 
i=1  i=2  i=1  i=2 


m 


m 


X  aii+2Xti-itiai.ij 


i=1 


i=2 


and 


(t-r-'y)  = 


m  m 

X  yi^^ii  +  X  (yi-i+yi)<Ji-i.i) 


i=1  i=2 


m  m 


X  tiyi«yii  +  X  (yiti-l+tjyi_i)aj_i,i) 

J=1  i=2 


8.  GRAPHICAL  ANALYSIS  OF  INTERVAL-CENSORED  DATA 


37 


Graphical  methods  are  used  to  identify  an  appropriate  distribution  to  model 

the  distribution  of  a  lifetime  variable  T.  Suppose  ti,  t2 .  tn  are  n  independent 

realizations  of  T  and  no  censoring  occurred.  One  can  plot  the  empirical 
distribution  function  of  the  data  on  probability  papers  such  as  Exponential, 
Gamma,  Weibull,  Lognormal,  and  Extreme  value  distribution  probability  paper.  If 
the  graph  of  the  empirical  distribution  function  of  the  data  resembles  a  straight  line 
on  a  particular  probability  paper,  the  corresponding  probability  model  is  used  to 
model  the  distribution  of  T. 

If  the  data  are  right-censored,  the  Kaplan-Meier  estimate  of  the  distribution 
function  of  T,  which  is  an  analogue  of  the  empirical  distribution,  is  used.  The 
Kaplan-Meier  estimate  can  be  plotted  on  a  variety  of  probability  papers  to  identify 
the  appropriate  parametric  distribution  to  model  the  distribution  of  T. 

For  interval-censored  data,  one  can  also  employ  probability  papers  to 
identify  the  underlying  parametric  model  for  the  lifetime  distribution.  Nelson  (1982) 
has  demonstrated  the  use  of  the  Weibull  probability  paper  for  plotting  interval- 
censored  data.  He  indicated  that  Exponential  probability  paper.  Lognormal 
probability  paper,  and  Extreme  value  probability  paper  can  be  used  analogously  to 
plot  interval-censored  data.  To  implement  the  procedure  of  Nelson,  one  needs  to 
create  all  these  probability  papers  before  plotting  the  data. 

In  many  statistical  software  packages  such  as  STATPAC  (Nelson  et  al., 
1978)  and  SAS  (1994),  creation  of  these  special  probability  papers  is  inherent  as 
a  part  of  model  building  endeavor.  In  this  section,  we  want  to  demonstrate  how 
regular  graph  paper  can  be  used  to  plot  interval-censored  data.  We  bring  Logistic 
and  Pareto  distributions  into  our  ambit  in  addition  to  the  distributions  considered 
by  Nelson  (1982). 

We  now  spell  out  the  problem  more  explicitly.  Suppose  we  want  to 
estimate  the  lifetime  distribution  T  of  a  product.  The  sampling  scheme  that  is  in 
focus  in  this  section  runs  as  follows.  Select  a  sample  of  n  units  of  the  product  and 

set  them  to  work.  Inspect  the  units  at  times  ti,  tz . U,  where  ti  <  ta  <  ...  <  U  are 

fixed.  Let  Xi  =  #  units  failed  during  the  interval  (tj-i,  tj],  i=1,  2 . m  and  Xm+i  =  # 

units  failed  during  (U,  «»),  where  (tj-i,  ti]  =  (0,  ti]  for  i  =  1.  Using  the  data  Xi,  Xz, 
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Xm+1,  appropriate  graphical  methods  will  be  developed  to  identify  the 
underlying  parametric  lifetime  distribution.  The  parametric  distributions 
considered  are  Exponential,  Weibull,  Lognormal,  Pareto,  Logistic  and  Extreme 
value  distribution. 


8.1  Graphical  Analysis  for  the  Exponential  Distribution 

The  probability  density  function  of  Exponential  distribution  is  given  by 


fT(t)  = 


ee 


-et 


t>0, 

otherwise. 


where  0  >  0  is  the  unknown  parameter  of  the  model.  Let 

Pi=P(ti_-,<T<ti), 

=  e-eti-i  i  =  i,  2,  ...,  m, 

and 

Pm+1  “  ^  ^m) 


=  e 


Note  that  for  any  j,  j  =  1 , 2, ....  m. 


ipi  =  1-e 

i=1 


-eti 


Let 


and 


Define 


m+1 

qj=  iPi. 

i=j+1 


rj  =^n(qj) 
m+1 

=  ^n(  XPi) 

i=j+l 

=  -0tj,  j  =  1,2,  ...,  m. 


m+1  Y 

Yj  =  ^n(  I  ^). 

i=j+1 
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The  joint  distribution  of  Xi,  X2,  Xm.,  Xm+i  is  multinomial  (n,  pi,  pa,  ....  pm,  Pm+i). 
.  m+1  X* 

Let  Xj  =  Note  that  Yj  =  ^n( — L).  Observe  that  X.-  is  distributed  as 

i=j+l 

Binomial  (n,  Pj) .  Using  the  A-method  (Rao,  1973,  p.  426),  we  get 

* 

E(Yj)  =  E(«n(^)) 
n 

•I  m+1 

s  ^n(-n  XPi) 

"  i=j+1 

=  -etj. 

which  can  be  written  as 

E(Yj)  =  ptj, 

where  p  =  -0  ,  ]  =  1,  2 . m.  It  is  now  clear  that  the  random  variables  Yi,  Y2 . 

Ym  conform  to  a  linear  model  through  the  origin.  This  observation  can  be 

exploited  as  follows.  Let  Xi,  xa,  ...,  Xm,  Xm+i  be  the  realizations  of  Xi,  X2 . Xm, 

m+1 

Xm+1,  respectively.  Let  y|  =  ^n(  ^  j  =1,2,  ...,  m.  Plot  (tj,  yj),  ]  =  1, 2 . m 

i=j+l  ^ 

on  a  regular  graph  paper.  If  the  points  lie  more  or  less  on  a  straight  line. 
Exponential  distribution  for  T  is  worth  a  try.  If  this  is  the  case,  the  preliminary 
estimate  of  0  can  be  obtained  from  the  graph.  This  estimate  is  given  by  the 
negative  slope  of  the  straight  line. 


8.2  Graphical  Analysis  for  the  Weibull  Distribution 

The  probability  density  function  of  Weibull  distribution  is  given  by 


0 


fT(t)  = 

where  A,  >  0  and  y  >  0  are  the  unknown  parameters  of  the  model.  Let 


t>0, 

otherwise. 


Pi=P(tj_-,<T<ti), 

=  i  =  1,2,  ....  m. 
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and 


Pm+1  “  >  ^m) 


=  e 


Note  that  for  any  j,  j  =  1 , 2 . m, 


Let 


and 


Define 


i 

SPi  =  1-e"*i. 

i=1 


m+1 

qj=  iPi. 

i=j+i 


fj  =  £n(-£n(qj) 
m+1 

=  £n(-£n  XPi) 
i=j+l 

=  ^nX,  +  YMj,  j  =  1,  2 . m. 


m+1  y 

Y|  =  ^n(-£n  X  — ). 


i=j+1 


n 


The  joint  distribution  of  Xi,  X2 . Xm,.  Xm+i  is  multinomial  (n,  pi,  pa,  ....  Pm,  Pm+i). 

.  m+1  X* 

Let  Xj  =  X^i-  Note  that  Yj  =  £n(-^n — ^).  Observe  that  X,  is  distributed  as 
i=j+i  n 

Binomial  (n,  qj).  Using  the  A-method,  we  get 

xJ 

E(Yj)  =  E(fnHn— !-) 

’  n 

•|  m+1 

=  ^n(-£n  -n  XPi) 
i=j+1 

=  in  A,  +  YMj, 

which  can  be  written  as 

E(Yj)  =  a  +  psj, 
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where  a  =  ^nX,,  p  =  y,  and  Sj  =^ntj ,  j  =  1 ,  2,  It  is  now  clear  that  the  random 
variables  Yi,  Ya,  Ym  conform  to  a  linear  model.  This  observation  can  be 
exploited  as  follows.  Let  Xi,  xa,  ....  Xm,  Xm+i  be  the  realizations  of  Xi,  Xa,  ....  Xm, 

m+1 

Xm+1,  respectively.  Let  yj  =  £nHn  ^  — ).  j  =1.2 . m.  Plot  (Mj,  yj),  j  =  1, 

i=j+1  ^ 

2, ...,  m  on  a  regular  graph  paper.  If  the  points  lie  more  or  less  on  a  straight  line, 
Weibull  distribution  for  T  is  worth  a  try.  If  this  is  the  case,  preliminary  estimates  of 
X  and  y  can  be  obtained  from  the  graph.  The  shape  parameter  y  is  estimated  by 
the  slope  of  the  straight  line,  and  the  scale  parameter  X  is  estimated  by  the 
exponent  of  the  intercept  of  the  straight  line. 

At  this  juncture,  it  is  appropriate  to  spell  out  the  essential  difference 
between  our  approach  and  Nelson’s  approach  in  plotting  interval-censored  data. 
We  and  Nelson  (1982)  use  the  same  relation 

^n(-^nP(T  >  tj))  =  ^nX  +  y^ntj,  j  =  1,  2,  ...,  m. 

Nelson  (1982)  takes  the  x-axis  on  the  logarithmic  scale,  y-axis  on  the  ^n(-^n) 

scale,  and  the  interval-censored  data  in  the  form  (tj,qj),  j  =  1,  2 . m  is  plotted 

m+1 

on  this  specially  created  paper,  where  q,  =  Y  — .  In  our  approach,  we  advocate 

i=j+i" 

plotting  (£ntj,£n(-^n  qj))  directly  on  a  standard  graph  paper.  There  are  two 

distinct  advantages  our  method  of  plotting  has  over  the  method  of  Nelson. 

1 .  The  data  points  (^ntj,^n(-£n  qj ),  j  =  1,  2,  ...,  m,can  be  plotted  precisely  on  the 

standard  graph  paper.  On  the  Weibull  probability  paper,  the  point  (tj,  qj)can 
only  be  identified  approximately. 

2.  Preliminary  estimates  of  the  scale  and  shape  parameters  of  the  Weibull 
distribution  are  easier  to  obtain  from  the  standard  graph  paper. 

As  an  illustration,  for  the  following  data  (Table  8.1)  on  167  identical  parts  in 
a  machine  (Nelson,  1982,  p.  415),  the  two  plots  are  presented  (  Figures  8.1  and 
8.2). 
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At  certain  ages,  the  parts  were  inspected  to  determine  which  had  cracked 
since  the  last  previous  inspection.  The  inspection  times  are  6.12,  19.92,  29.64, 
35.40,  39.72,  45.24,  52.32,  and  63.48  months  respectively.  The  data  appear  in 
Table  8.1;  it  shows  the  months,  in  service  at  the  start  and  end  of  each  inspection 
period  and  the  number  of  cracked  parts  found  in  each  period.  For  example, 
between  6.12  and  19.92  months,  16  parts  cracked.  Seventy  three  parts  survived 
the  latest  inspection  at  63.48  months. 

Table  8.1 

Part  Cracking  Data 


Inspection  Intervals 

Number 

Cracked  Cumulative 

■ 

q j 

(0.00, 

6.121 

5 

5 

0.970 

0.023 

-3.493 

(6.12, 

19.921 

16 

21 

0.874 

0.126 

-2.007 

(19.92, 

29.641 

12 

33 

0.802 

0.198 

-1.513 

(29.64, 

35.401 

18 

51 

0.695 

0.305 

-1.009 

(35.40, 

39.721 

18 

69 

0.587 

0.413 

-0.629 

(39.72, 

45.241 

2 

71 

0.575 

0.425 

-0.591 

(45.24, 

52.321 

6 

77 

0.539 

0.461 

-0.481 

(52.32, 

63.481 

17 

94 

0.437 

0.563 

-0.189 

(63.48, 

_ _ 

73 

167 

Total 

167 
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Figure  8.1 

Estimation  of  the  Weibull  Distribution  Parameters  for  the 
Part  Cracking  Data  Using  Standard  Graph  Paper 

The  preliminary  estimates  of  the  Weibull  parameters  from  the  standard 
graph  paper  are  i  =  0.00452  and  y  =  1 .636. 
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Figure  8.2 

Estimation  of  the  Weibull  Distribution  Parameters  for  the 
Part  Cracking  Data  Using  Weibull  Probability  paper 


The  preliminary  estimates  of  the  Weibull  parameters  from  the  Weibull 
Probability  paper  are  i  =  0.0017  and  y  =  1 .49. 


45 


8.3  Graphical  Analysis  for  the  Lognormal  Distribution 

The  probability  density  function  of  the  Lognormal  distribution  is  given  by 


fT(t)  = 


1 


taV^ 

0 


20*^ 


t>0, 

Otherwise, 


where  a  >  0  and  -o®  <  p.  <  oo  are  the  unknown  parameters  of  the  model.  Let 

Pi=P(ti_i<T<ti), 


and 


where 


=  <& 


-<D 

\  J  [ 


,  i  =  1,  2,  ....  m, 


Pm+1  ~  P(T  >  tpn) 


<I>(x)  =  I 


^  e  2  dz. 


Note  that  for  any  j,  j  =  1 , 2, ....  m, 


I 


Let 


and 


SPi  = 

i=1 


f^ntj-p^ 


V  ^  y 


qj  =  XPi  =  1-<I>  ^ 

i=j+i 


V  ^  y 


rj  =0)  \l-qj), 

=  £ntj,  i  =  1,  2,  m. 

cr  o  * 


Define 
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The  joint  distribution  of  Xi,  X2,  Xm.,  X^+i  is  multinomial  (n,  pi,  pz,  p^,  pm+i). 

‘  "1+1  xj 

Let  Xj  =  Note  that  Yj  =4)  '(1 — -).  Observe  that  Xj  is  distributed  as 

i=j+l  ^ 

Binomial  (n,  qj) .  Using  the  A-method,  we  get 


which  can  be  written  as 


E(Yi)  =  E4>-'(1-^i-) 
n 

1  "1+1 

=  ^  IPi) 

”  i=j+1 
1  -  X 

=  — Mi, 

a  a  ' 


E(Yj)  =  a  +  pSj, 

where  a  =-p/a,  p  =  1/a,  and  Sj  =^ntj ,  j  =  1,  2,  ...,m.  It  is  now  clear  that  the 

random  variables  Yi,  Y2,  ....  Ym  conform  to  a  linear  model.  Let  Xi,  X2 . Xm,  Xm+i 

be  the  realizations  of  Xi,  X2,  ....  Xm,  Xm+i,  respectively.  Let 

_1  *^+1 X 

yj  =0  '(1-  X  — ).  j  =1.2,  ...,  m.  Plot  (Mj,  yj),  j  =  1,  2,  ...,  m  on  a  regular 
i=j+i  1 

graph  paper.  If  the  points  lie  more  or  less  on  a  straight  line.  Lognormal 
distribution  for  T  is  a  possible  model.  If  this  is  the  case,  preliminary  estimates  of  jx 
and  a  can  be  obtained  from  the  graph.  The  estimate  a  of  a  is  given  by  1/p,  and 
the  estimate  ji  of  p  is  given  by  -o/p,  where  a  and  p  are  the  intercept  and  slope, 
respectively,  of  the  free  hand  straight  line  passing  close  to  the  points  plotted. 


8.4  Graphical  Analysis  for  the  Pareto  Distribution 

The  probability  density  function  of  Pareto  distribution  is  given  by 


fT(t)  = 


yaV 

^  tY+1 
0 


t  >a, 

otherwise, 


where  y  >  0  is  the  unknown  parameter  of  the  model.  Assume  that  a  is  known  In 
this  case,  the  first  inspection  interval  is  (a,  ti] .  Let 

Pi  ~  "I"  —  fi)i 


and 


Pm+1  ~  ^  ^m) 


Vm  } 


Note  that  for  any  j,  j  =  1 , 2, ...,  m. 


Let 


and 


Define 


ipi  =  1- 

i=1 


k'i; 


m+1 

qj=  iPi 

i=j+1 


v'i; 


m+1 

rj  =^nqj  =^n(  XPi). 

i=j+1 

=  j  =  1.  2 . m. 

*j 
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The  joint  distribution  of  Xi,  X2,  X^.,  Xm+i  is  multinomial  (n,  pi,  pa,  ....  pm,  Pm+i)- 

Let  Xj  =  Note  that  Yj  =£n—.  Observe  that  Xj  is  distributed  as 

i=]+l  1 

Binomial  (n,  qj) .  Using  the  A-method,  we  get 

xJ 

E(Yj)  =  E(«n-L) 
n 

-•  m+1 

=  ^n(-n  XPi) 
i=j+i 

=  Y^n(-^),  1  =  1,2 . m. 

It  is  now  clear  that  the  random  variables  Yi,  Y2,  ....  Ym  conform  to  a  linear  model 
through  the  origin..  Let  Xi,  X2,  ....  Xm,  Xm+i  be  the  realizations  of  Xi,  X2 . Xm, 

m  +  1  g 

Xm+1,  respectively.  Let  yj  =  £n(  X  — ).  j  =1.  2,  ....  m.  Plot  (^(-).  yj),  j  =  1,  2, 

i=j+i  'I  1] 

....  m  on  a  regular  graph  paper.  If  the  points  lie  more  or  less  on  a  straight  line, 
Pareto  distribution  is  a  possible  model  for  T.  If  this  is  the  case,  preliminary 
estimate  y  of  y  can  be  obtained  from  the  graph.  This  estimate  is  given  by  the 
slope  of  the  straight  line. 


8.5  Graphical  Analysis  for  the  Logistic  Distribution 

The  probability  density  function  of  Logistic  distribution  is  given  by 

'  XytY-i 


fT(t)  = 


O  +  XX'^f 
0 


t>0, 

otherwise. 


where  X>  0  and  y  >  0  are  the  unknown  parameters  of  the  model.  Let 

Pi  =P(ti_i  <T<ti), 


1 


1 


1  +  X.tj^.j  1  +  X,tj^ 


,  i  =  1,  2,  ...,  m. 


and 
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Let 


and 


Pm+1  “  >  ^m) 

1 


Note  that  for  any  j,  j  =  1 , 2 . m, 


Define 


1  +  A.t^, 


1 


i=1  1  +  X,t  i 


m+1  -I 

qi=  Spi  = 


i4;i  i+wr 


1-qj 

r,  =  £n  ( - i-)  =  ^n 

qj  ' 


=  ^n  +  Y^n  tj,  j  =  1,  2 . m. 


^  1 
Ipi 

i=1 

=  £n 

f  \ 

1  1 

m+1 

m+1 

IPi 

V'=j+1  J 

i^i=j+1  J 

Yj  =  £n| 


1 


m+1  Y 

I  ^ 


-1 


Li=j+1 


n 


The  joint  distribution  of  Xi,  Xa . Xm.,  Xm+i  is  multinomial  (n,  pi,  pa . Pm,  Pm+i). 


m+1 


Let  Xj  =  ]^Xj.  Note  that  Yj  =^n(-!V-1)-  Observe  that  Xj  is  distributed  as 

i=j+i  Xj 


Binomial  (n,  qj) .  Using  the  A-method,  we  get 
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E(Yj)  =  E£n(-L-1) 

Xj 


m+1 

n  IPi 

i=j+i 
=  ^nX  +  yMj, 

which  can  be  written  as 


E(Yj)  =  a  +  pSj, 

where  a  =  (nX,  p  =  y  and  Sj  =^ntj ,  j  =  1,  2,  It  is  now  clear  that  the  random 

variables  Yi,  Ya . Ym  conform  to  a  linear  model.  We  use  this  as  follows.  Let  Xi, 

X2,  ....  Xm+1  be  the  realizations  of  Xi,  Xa,  ....  Xm,  Xm+i,  respectively.  Let 

-1 

graph  paper.  If  the  points  lie  more  or  less  on  a  straight  line,  Logistic  distribution 
for  T  is  worth  a  try.  If  this  is  the  case,  preliminary  estimates  of  y  and  a  can  be 

obtained  from  the  graph.  The  estimate  y  of  y  is  given  by  p,  and  the  estimate  X  of 

X  is  given  by  e“,  where  a  and  p  are  the  intercept  and  the  slope,  respectively,  of 
the  free  hand  straight  line  passing  close  to  the  points  plotted. 


8.6  Graphical  Analysis  for  the  Extreme  Value  Distribution 

The  probability  density  function  of  Extreme  value  distribution  is  given  by 

^x  (x)  =  y  -  oo  <  X  <  oo 

where  -oo  <  e  <  oo  and  y  >  0  are  the  unknown  parameters  of  the  model.  Let 

T  =  e^. 

It  can  be  shown  that  the  probability  density  function  of  T  is  given  by 
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where  A,  =  e“^,  which  is  the  probability  density  function  of  the  Weibull 
distribution  with  parameters  y  and  X,.  Therefore,  to  estimate  the  parameters  of  the 

Extreme  value  distribution,  we  fit  the  Weibull  distribution  to  the  T=e^  variable. 
The  estimate y  of  y  is  equal  toy, and  the  estimate  9  of  0  is  equal  to  -fnX/y . 

9.  LINEAR  MODEL  APPROACH  TO  INTERVAL-CENSORED  DATA 

Suppose  we  want  to  estimate  the  lifetime  distribution  T  of  a  product.  The 
sampling  scheme  that  is  the  general  periodic  inspection  plan  outlined  in  Section  8. 
Using  the  data  Xi,  X2,  ...,  Xm,  Xm+i,  the  main  focus  in  this  section  is  on  the 
estimation  of  the  parameters  of  the  underlying  distribution  of  T  following  the  Linear 
Model  Approach.  This  method  of  estimation  is  simpler  than  the  method  of 
maximum  likelihood.  These  two  methods  give  asymptotically  consistent 
estimators.  We  examine  how  they  perform  in  small  samples.  The  parametric 
distributions  considered  are  Exponential,  Weibull,  Lognormal,  Pareto,  Logistic, 
and  Extreme  value  distributions. 

9.1  Linear  Model  Approach  for  the  Exponential  Distribution 

The  starting  point  is  the  collection  of  random  variables, 

m+1  V. 

Yj  =£n(  — ).  i  =  1,  2 . m. 

Using  the  A-method,  we  have  noted  that 

x[ 

E(Y|)  =  E«n(-L)) 
n 

s  ^n(-nqj) 
n  * 

=  -etj, 

which  can  be  written  as 

E(Y,)  =  pt,. 
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where  p  =  -0  ,  j  =  1 , 2,  Using  the  A-method  again,  it  can  be  shown  that 


Var(Y|)  =  Var(^n(— L)) 


n 


and  for  j  <  k, 


Cov(Yj.  Yk)  =  Cov(£n(-^),^n(^)) 

n  n 


=  Var(Xj) 

nqj 


.  n2 

d  Xi 

-^^n(-L) 

dX=  n 


XpEXi 


=  Cov(X|.Xk) 


- J 

Q. 

_ » 

[dXj  n  J 

L«Xk  "  J 

Xk=EXi, 

x;=EXj 


m+1  m+1 

=  Cov(  XXj,  IX|,)| 

i=j+1  i=k+1 


d  Xj 

.  ^n( 

dXj  n 

[dXk  "  J 

Xi=EX, 


m+1 

k 

m+1 

H  Xi 

Var(  2X0  + 

I 

SCov(Xr.Xs) 

i=k+1 

r=j+1s 

=k+1 

dXj  n 

d  .  /Xk, 


dXi 


Xk=EXj, 

X;=EXj 


k  m+1 

nqk(1-qk)+  S  Z-nPrPs 

r=j+1s=k+1 


1  1 


nqj  npk 


1- 


dj 


nqj 

In  general,  for  j  ^  k, 


Cov(Yi,  Yk)  =  ^ 
nqs 


where  s  =  Min  (ji  k). 

The  variances  and  covariances  of  Y|’s  have  the  following  properties: 

Var(Yi)  =  Cov(Yi.  Ya)  =  Cov(Yi,  Y3)  =  ...  =  Cov(Yi,  Ym) : 
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Var(Y2)  =  Cov(Y2,  Y3)  =  Cov(Y2,  Y4)  =  ...  =  Cov(Y2,  Y^,) ; 


Var(Y,T,.i)=Cov(Y,n.i,Ym). 

Let  us  introduce  a  special  notation  exploiting  the  special  structure  of  these 
variances  and  covariances.  Let 


e|<=Var(Yk)  =  -H!i-,  k  =  1,2,  ....  m. 

nqk 

These  computations  can  be  summarized  as  follows: 


E(Y)  =  E( 


fYi^ 

r-M 

Y2 

C\J 

I 

Y3 

)  = 

I 

CO 

• 

• 

e  =  tp, 


the  dispersion  matrix  Z  of  Y  is  given  by 

Gf  0.J  0^  0.| 

02  02  02 
03  03  03 

04  04 

E  =  I  ®5 


01 

02 

03 

04 

05 


01 

02 

03 

04 

05 


0(m-1)  0(m-1) 

0m  j 


Let  Z  ^  =  (ojj).  In  Section  6,  we  have  identify  the  inverse  Z  ^  to  be  a  tri-diagonal 
matrix.  If  Z  is  known,  the  best  linear  unbiased  estimator  Gp  of  0  is  given  by 
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0n 


m  m 

I  +  I  (tMyi+tiyi_i)CTi_-,j 

_ 152 _ 


S  tfoii  + 

m 

i=1 

i=2 

Since  Z  is  unknown,  one  can  estimate  Z  consistently  by  Z 

,  where 

x  A  A 

e,  _e, 

§1  §1  01 

...  01 

ei  ' 

§2 

§2  §2  §2 

...  §2 

02 

§3  03  §3 

...  §3 

03 

§4  §4 

...  04 

04 

Z  = 

...  §5 

05 

1 

A 

0(m-1) 

A_ 

and 


The  modified  estimator  0n  of  0  is  given  by 

0n  =(t'i:''’t)"hl"‘'Y 
m  m 

I  +  X  (tMyi+tiyi-i)Vi.i 

_i=i _ £2 _ 

m  m 

X  tfaii+2£  (tMti)aj_ij 
i=1  i=2 

where 

i-’  =  (Oij). 

Some  useful  facts  emerge  from  these  deliberations. 
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1 .  We  have  an  explicit  formula  for  the  estimator  of  0.  The  maximum  likelihood 
method  does  not  provide  an  explicit  formula  for  the  maximum  likelihood 
estimator  of  0. 

2.  The  new  estimator  is  also  consistent.  For 

plim  0n  =  plim(t'i‘‘'t)"‘'t'i"‘'Y 
n— 

=  (t'  plim  i’^t)"^t'(plim  t"^)(plim  Y),  (Slutsky' s  theorem) 

n— >oo  n— n— 

=  (t'2:'^t)“''t'z"^tp 
=  p  =  0. 

In  one  special  case,  the  maximum  likelihood  method  provides  an  explicit 
formula  for  its  estimator.  This  estimator  differs  from  the  linear  model  estimator. 
Let  t|  =  ito,  i  =1 , 2, ...,  m  for  some  to  >  0. 

Linear  model  estimator: 


m  m 


0n  = 


i=1 


i=2 


m  m 


i=1 


Maximum  likelihood  estimator: 


0n  = — £n| 


to 


i=2 


f  m  ^ 

mn+  Z(i-nri)Xj 

i=i _ 

m 

mn+  2(i-1-m)Xj 
i=1 


In  the  particular  case  m  =  2,  one  can  write  down  the  estimators  in  terms  of  the 
data  Xi,  X2,  and  X3. 

Linear  model  estimator: 


0n  =  T^( 


XiX 


/X2  +  X3  . 


nX. 


to  (Xi+X2)(X2+X3)  '  X3  '  (Xi+X2)(X2  +  X3) 

Maximum  likelihood  estimator: 


.X2  +  X3 .. 


n 


1 

0n=— ^n 


2n-X 


1 


to  2n-2Xi-X2 
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A  warning  should  be  sounded  at  this  juncture.  In  practice,  especially  if  n  is 
small,  the  observed  values  of  some  Xj’s  could  be  zero.  In  this  case,  the  estimated 

dispersion  matrix  i  is  singular.  The  linear  model  estimator  does  not  exist.  The 
maximum  likelihood  is  also  not  free  of  problems.  The  iterative  procedure  used  to 
estimate  the  parameters  fails  to  converge.  One  can  note  that  P(Xi  =  0)  converges 
to  0  exponentially  fast  as  n  -^  <».  Consequently,  the  occurrence  of  the  event  { Xi  = 
0}  is  unlikely. 


9.2  Linear  Model  Approach  for  the  Weibull  Distribution 

We  begin  with  the  random  variables, 


m+1  y 

Yj  =  ^n(-£n  I  ^). 
i=j+i 


Using  the  A-method,  we  have  noted  that 

E(Yj)  =  ^n  X  +  Y^ntj, 

which  can  be  written  as 


E(Yj)  =  a  +  psj, 


where  a  =  A,,  p  =  y,  and  Sj  =^nt] ,  j  =  1 ,  2,  ...,m.  It  is  now  clear  that  the  random 

variables  Yi,  Y2,  ...,  Ym  conform  to  a  linear  model.  Using  the  A-method  again,  it 
can  be  shown  that 


Var(Yj)  =  Var(^n(-^n(— )) 


n 


=  Var(Xj) 


dXi 


I 

-£n(-^n(— 


n 


Xi=EX, 


nqj(^nqj) 


2  ’ 


and  for  j  <  k. 
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Xi 


Cov(Yj,  Yk)  =  Cov(^n(-^n(^),^n(-^n(^)) 

n  n 


=  Cov(Xj.Xk) 
1-qj 


d  ,  X| 

d  xl 

dXj  n 

[dXk  "  J 

Xk=EXk 

x;=Ex; 


nqj(^nqj)(^nqk) 
In  general,  for  j  it  k, 


Cov(Yi,  Yk)  = - - , 

'  nqs(£nqj)(^nqk) 


where  s  =  Min  (j,  k). 

Let  us  introduce  a  special  notation  exploiting  the  special  structure  of  these 
variances  and  covariances.  Let 

1-qj 


9]  = 


nqj 


and 


7Cj  = 


^nqi 


Then, 

and  for  j  it  k. 


Var(Yj)  =  e|Kj,  j  =  1,  2,  ...,  m, 
Cov(Yj,  Yk )  =  SgTtjTtki 


where  s  =  Min  (j,  k). 

These  computations  can  be  summarized  as  follows: 


fYO 

er^x■^  ^ 

Y2 

1  £n\2 

E(Y)  =  E( 

Y3 

)  = 

1  ^nt3 

'.Ym; 

^1  ^ntfpy 

=  sp. 


58 


the  dispersion  matrix  Z  of  Y  is  given  by 


71^01  7IiK201  7ti7l30i  3ri7l40i 


z  = 


^202  7t2^3®2  7l27t402 

2 

^4®4 


^2^(m-1)®2  ^2^m®2 

^3^(m-1)®3  ^3^01^3 


2 

^(m-1)®{fn-1)  ^(m-1)^m®(m-1) 


^m®m 


Let  Z  =  (cTij).  The  inverse  Z  Ms  a  tri-diagonal  matrix.  If  Z  is  known,  the  best 
linear  unbiased  estimators  a  and  p  of  a  and  p,  respectively,  are  given  by 

=  (S'Z“‘'S)~^S'Z“''y. 


.-1 


Since  Z  is  unknown,  one  can  estimate  Z  consistently  by  Z ,  where 


Z  = 


^2A 


71^01  7IiK201  7Ci7t30i  7Ci7l40i 


7t|02  7t2^3®2  • 

71^03  71371403  . 

^464  . 


^2^(m-1)®2 

^3^(171-1)03 


^2^m02 

*3^'m03 


^(m-i)0(m-1)  ^(m-1)^m0(m-1) 


0j  = 


n-xi 


nxi 


and 


7Cj  = 


in- 


n 
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The  modified  estimators  of  d  and  p  of  a  and  p  are  given  by 

d"\  .  .  1-1 

pj 

Ja  bVVdi 

'.b  cj  l^d2 

where 

m  m 

a  =  X  aii  +  2Xaj_i.j, 

i=1  i=2 

m  m 

b  =  Xtidii  +  ]£  (tj-i 
i=1  i=2 

m  m 

i=1  i=2 

m  m 

di  =  S  Vian  +  X  (yM  +  yi)ai_i,i, 
i=1  i=2 

m  m 

d2  =  Z  tiyi^ii  +  X  (yjti-i  +  tiyi_i)ai_ij. 
i=1  i=2 

Therefore, 

.  cdi-bdo 
a  =  — ' - 

ac-b^ 

and 

-  _  ad2  -bdi 
ac-b^ 

As  an  illustration,  we  fit  the  Weibull  distribution  to  the  Part  Cracking  Data 
given  in  Table  8.1.  Table  9.1  presents  the  observed  cracked  parts  and  the 
expected  cracked  parts  under  the  maximum  likelihood  and  linear  model  methods. 
Under  the  linear  model  method,  the  distribution  of  the  lifetime  T  is  Weibull 
(0.0081,  1.1037):  and  under  the  maximum  likelihood  method,  the  distribution  of  T 
is  Weibull  (0.00178,  1.4825). 
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Table  9.1 


li 


Fitting  Weibull  Distribution  to  the  Part  Cracking  Data 


Inspection 

Intervals 

Obsen/ed 
Number  cracked 

GRiSiSI 

Expected 
Number  cracked 
(Maximum  Likelihood) 

OSI 

5 

9.71 

4.30 

16 

23.32 

18.96 

12 

15.31 

16.30 

Essn 

18 

8.46 

9.96 

18 

6.04 

7.44 

8 

16.09 

20.88 

17 

12.43 

16.77 

■m 

73 

75.64 

72.40 

Total 

167 

167 

167 

The  chi-square  goodness  of  fit  values  are 
=45.578  (p<  0.001), 
and 

=  31142  (p<  0.001). 

The  chi-squared  values  are  significant.  The  null  hypothesis  that  the  underlying 
distribution  is  Weibull  is  rejected. 


9.3  Linear  Model  Approach  for  the  Lognormal  Distribution 

Note  that  the  random  variables 

^  m+1  V. 

Yj=<I>"^1-  I  j  =  1,2,  ....  m. 
i=j+i  "" 


have  expectations, 


E(Yj)  =  E<I>-^(1 - i-) 

*  n 

=  -— +  -lMj, 
a  a  ' 


which  can  be  written  as 
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E(Yj)  =  a  +  pSj, 

where  a  =-|j/ct,  p  =  1/a,  and  Sj =^n  tj ,  j  =  1 , 2,  Using  the  A-method  again,  we 
get 


Xi 


Var(Yj)  =  Var(a)"^(1 - i-)) 

^  n 


=  Var(Xj) 


dX=  n 


.  n2 


‘J 


X,=EXj 


Thus, 

.271.,  .  k'^t-qj)]^ 

Var(Yj)s— qj(1-qj)eL  'J, 

and  for  ]  <  k. 


Cov(Yj,  Yk)  =  Cov 

i).  *-'(1-4) 

=cov(x;,x;) 

‘'.‘■•-'(I-’'')! 
dXj  n  J 

[dXk  "  J 

X,=EX, 

[<l>-\l-q))f  [o-yi-qK)]^ 

=  -^qk(1-qj)e  2  e  2  , 

In  general,  for  j  vt  k, 

[‘l>-\l-qj)f  [‘I>''(1-qi<)f 

Cov(Y|.  Yk)  =  .^qs(1-q,)e'  2  e“  2  , 

where  r  =  Min  (j,  k)  and  s  =  Max  (j,  k). 

Let  us  introduce  a  special  notation  exploiting  the  special  structure  of  these 
variances  and  covariances.  Let 


62 


2k 


and 


ej=^(1-qj)e  2 


~  ^k®  ^ 


Then, 


and  for  j  9^  k, 


Var(Yj)  =  0j7tj,  j  =  1,  2,  ....  m, 


Cov(Yj,Yk)  =  er7ts, 
where  r  =  Min  (j,  k)  and  s  =  Max  G,  k). 

These  computations  can  be  summarized  as  follows: 


E(Y)  =  E( 


ri 

CM 

>- 

CM 

C 

Y3 

II 

1  M3 

a 


=  sp. 


the  dispersion  matrix  Z  of  Y  is  given  by 


Z  = 


7Ci0i  71201 

7C301 

7^401 

7t(nri-1)01 

71202 

71302 

K402 

•  •  7I(m_i)02 

7^77102 

71303 

71403 

•  •  7t(nri_i)03 

7177103 

7C404  . 

.  . 

• 

Using  the  data  one  can  estimate  the  dispersion  matrix, 
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z  = 


fti0i  71:201 

71301 

71401  •  •  • 

71202 

CM 

CO  • 

CO 

7t402  .  .  . 

*(m-1)62 

7C303 

71403  •  •  ■ 
71:404  ,  .  . 

^(m-1)03 

^m03 

A  /y  ^ 


^m®m 


where 


and 


2k 


•  i2 


Xj  - 


0.  =  -(1_:i)e  2 
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n 


\  0 
n 

If  Z  is  known,  the  Best  Linear  Unbiased  Estimators  dand  pof  a  and  p  are  given 
by 


P 


=  (S'Z-''S)"''SZ"''y 


Since  Z  is  unknown,  the  estimate  Z  of  Z  is  used,  i.e., 

=  (sz-''s)-‘'sz-''y. 


^d^ 


vPy 


94  Linear  Model  Approach  for  the  Pareto  Distribution 

Note  that  the  random  variables 
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have  expectations, 


i=j+1  ” 


E(Yj)  =  ^n(f). 


which  can  be  written  as 

E(Vi)  =  PSj. 

where  p  =  y,  Sj  =  £n(— ) ,  j  =  1 , 2,  Using  the  A-method  again,  we 

I  i 


Var(Yj)  =  Var(^n(-i-)) 


.  i2 

d  Xi 


dXj  n 


XfEx; 


and  for  j  <  k, 


=  Var(X;) 
nqj  ’ 


X*  y* 

Cov(Yj,  Yk )  =  Cov(£n(^).  £n(-^)) 

n  n 


nqj  ■ 

In  general,  for  j  k, 

Cov{Yj,  Y„)  =  ^, 
nqs 

where  s  =  Min  (j,  k). 

The  variances  and  covariances  of  Yi’s  have  the  following  properties: 

Var(Yi)  =  Cov(Yi,  Y2)  =  Cov(Yi,  Y3)  =  ...  =  Cov(Yi,  Ym) ; 
Var(Y2)  =  Cov(Y2,  Y3)  =  Cov(Y2,  Y4)  =  ...  =  Cov(Y2,  Ym) : 


Var(Y,T,.i)  =  Cov(Ym.i,Y„,). 


Let  us  introduce  a  special  notation  exploiting  the  special  structure  of  these 
variances  and  covariances.  Let 


0k=Var(Yk)  =  -i^.  I<  =  1.2,  ...  m. 

nqk 

These  computations  can  be  summarized  as  follows: 


the  dispersion  matrix  I  of  Y  is  given  by 

^0-(  0-j  0-|  0-|  0-| 

02  02 

03  03  03 

04  04 


Using  the  data,  one  can  estimate  the  dispersion  matrix, 
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where 


0j  = 


n-x 


nxi 


Let  S  ^  =  (ay).  The  inverse  Z  ^  is  a  tri-diagonal  matrix.  If  Z  is  known,  the  best 
linear  unbiased  estimator  p  of  p  is  given  by 

P  =  (S'Z-''S)"‘'S'Z"''y 

Y  Yi^ii  I  Y  /  Yi  .  Yi-1  v 
2.-^+  1  (t— 

_  i=i  h  j=2  ‘i-1  M 


m 


m 


I  ^+2i; 


i=1  i=2 

Since  Z  is  unknown,  one  can  estimate  Z  consistently  by  Z ,  where 
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which  can  be  written  as 


E(Yj)  =  E<n(-!V-1) 

=  nA,  +  YMj, 


E{Yj)  =  a  +  psj, 

where  a  =  ^n  A,,  and  p  =  Y,and  Sj  =^ntj ,  j  =  1,  2,  Using  the  A-method  again, 
we  get 


Var(Yj)  =  Var(^n(— -1)) 


=  Var(X;) 
1 


Jx;=Ex; 


nqj(1-qj)’ 


and  for  j  <  k, 


Cov(  Yj ,  Y|, )  =  Cov(^n(-!V  - 1),  - 1))) 

X|  Xk 


=  Cov(x;,XO 


dX|  X| 


-^n(— -1) 


idXk  Xk 


Xk=EXk 


Xi=EX 


m+1  m+1 

=  Cov(  £X|,  2Xk) 

i=j+1  i=k+1 


-n 


Xj(n-Xj) 


-n 


IXk(n-Xk) 


Xk=EXk 

x;=Ex; 


k  m+1 

nqk(1-qk)+  X  X{-nPrPs) 

r=j+1s=k+1 

1 


T  -n  1 

-n 

1 _ 

jO 

I 

1 _ 

_nqk(n-nqk)_ 

nqj(1-qk) 

In  general,  for  i  k 


Cov(Yi,Yk)  = 


nqrO-<ls)’ 


where  r  =  Min  Q,  k)  and  s  =  Max  (j,  K). 
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Let  us  introduce  a  special  notation  exploiting  the  special  structure  of  these 
variances  and  covariances.  Let 


0j  = 


and 


nqj 


1 


Then, 


1-% 

Var(Yj)  =  0jKj,  j  =  1,  2,  ....  m, 


and  forj^tk, 

Cov(Y|,Yk)  =  0r7Cs, 
where  r  =  Min  (j,  k)  and  s  =  Max  O',  k). 

These  computations  can  be  summarized  as  follows: 


E(Y)  =  E( 


ri 

Y2 

CVJ 

c 

Y3 

)  = 

CO 

c 

• 

• 

^YmJ 

3 

a 


=  sp, 


the  dispersion  matrix  I  of  Y  is  given  by 

7tl0i 

7C202  ^402 

^3®3  ^4®3  • 

71404  • 

1 

2  = 


^(m-1)Q2 

^(m-1)03 


^tm02 
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Using  the  data  one  can  estimate  the  dispersion  matrix. 
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7ti0i  7t201 

7U301 
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7^202 
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^3^3 

71403  •  • 

^(m-1)03 

<CD 

.  . 

, 

where 


and 

n 

= - T- 

n-Xk 

If  Z  is  known,  the  Best  Linear  Unbiased  Estimators  d  and  pof  a  and  p  are  given 


by 


^d^ 


Since  Z  is  unknown,  the  estimate  t  of  Z  is  used,  i.e., 

=  (S'Z-''S')'''S'Z-''y. 


9.6  Linear  Model  Approach  for  the  Extreme  Value  Distribution 

The  probability  density  function  of  Extreme  value  distribution  is  given  by 

where  -oo  <  0  <  oo  and  y  >  0  are  the  unknown  parameters  of  the  model.  Let 
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T  =  e^. 


It  can  be  shown  that  the  probability  density  function  of  T  is  given  by 

f-r(t)  =  >.YtY*'’e'^*^'\  0<t<oo. 

where  X.  =  e~^,  which  is  the  probability  density  function  of  the  Weibull 
distribution  with  parameters  y  and  X.  Therefore,  to  find  the  linear  model 
estimators  of  the  Extreme  value  distribution  parameters,  we  fit  the  Weibull 

distribution  to  the  T=e^  variable.  The  estimate  y  of  y  is  equal  to  y,and  the 
estimate  0  of  0  is  equal  to  -£nX/y. 
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